1. 非線形の成長曲線モデル

• polynominal model
• latetn basis model
• piecewise models

2.多変量の成長曲線モデル

• parallel process growth model
• autoregressive latent trajectory Model
• latent change score model

２次(以上)の曲線で

polynominal model

データドリブンに

latent basis model

ある時点から傾きが異なる

piecewise model

ポジ感情の経時データ

Garland, E. L., Geschwind, N., Peeters, F., & Wichers, M. (2015). Mindfulness training promotes upward spirals of positive affect and cognition: multilevel and autoregressive latent trajectory modeling analyses. Frontiers in psychology, 6.

lavaanでまじめに書く

level =~ 1* bmi1 +1* bmi2 +1* bmi3 +
         1* bmi4 +1* bmi5 +1* bmi6 
slope =~ 0 * bmi1 + 1 * bmi2 + 2 * bmi3 + 
         3 * bmi4 + 4 * bmi5 + 5 * bmi6 
bmi1 ~~(vare)* bmi1 
bmi2 ~~(vare)* bmi2 
bmi3 ~~(vare)* bmi3 
bmi4 ~~(vare)* bmi4 
bmi5 ~~(vare)* bmi5 
bmi6 ~~(vare)* bmi6 

結構めんどい…

救世主登場

RAMpathパッケージ

• 潜在成長曲線系のlavaanコードを自動生成して実行してくれる関数が充実
  • 成長曲線モデル: ramLCM
    
  • 潜在差得点モデル:ramLCS
    
  • 2変量の潜在差得点モデル: ramBLCS
• lavaanのモデルを吐き出せるので、RAMpathで実行した後lavaanで微調整という流れで効率よくモデリング

ramLCM

• 切片のみのモデル (model=‘no’)

• 線形モデル (model=‘linear’)

• 二次曲線モデル (model=‘quadratic’)

• latent basisモデル (model = ‘latent’)

一行で簡潔に推定

一気にどん

library(RAMpath)
fit.all<-ramLCM(data=data,outcome=1:6, model='all')

個別にどん

fit.no<-ramLCM(data=data,outcome=1:6, model='no')
fit.linear<-ramLCM(data=data,outcome=1:6, model='linear')
fit.quadratic<-ramLCM(data=data,outcome=1:6, model='quadratic')
fit.latent<-ramLCM(data=data,outcome=1:6, model='latent')

切片のみモデル

lavaanコード

cat(fit.all$model$no)

## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6

• 傾き因子を排除
• 初期値から得点の変化なし

線形モデル

lavaanコード

cat(fit.all$model$linear)

## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  slope =~  0 * X1 + 1 * X2 + 2 * X3 + 3 * X4 + 4 * X5 + 5 * X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6

• 傾き因子の因子負荷が単調増加 (0, 1, 2, 3, 4, 5)

二次曲線モデル

lavaanコード

cat(fit.all$model$quadratic)

## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  slope =~  0 * X1 + 1 * X2 + 2 * X3 + 3 * X4 + 4 * X5 + 5 * X6 
##  quadratic =~  0 * X1 + 1 * X2 + 4 * X3 + 9 * X4 + 16 * X5 + 25 * X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6

• 二つ目の傾き因子の負荷は, 一つ目の傾き因子の二乗 (0, 1, 4, 9, 16, 25)

latent basisモデル

lavaanコード

cat(fit.all$model$latent)

## level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
##  slope =~  0 * X1 +start( 1 )* X2 +start( 2 )* X3 +start( 3 )* X4 +start( 4 )* X5 + 5 * X6 
##  X1 ~~(vare)* X1 
##  X2 ~~(vare)* X2 
##  X3 ~~(vare)* X3 
##  X4 ~~(vare)* X4 
##  X5 ~~(vare)* X5 
##  X6 ~~(vare)* X6

• 傾きの因子負荷を自由推定する
• データドリブンな遷移パターン -ramLCMの自動コードは若干ミスってる

コードミスを修正

modelXlat<-'
level =~ 1* X1 +1* X2 +1* X3 +1* X4 +1* X5 +1* X6 
slope =~  0 * X1 +start( 1 )* X2 +start( 2 )* X3 +start( 3 )* X4 +start( 4 )* X5 + 5 * X6 
X1 ~~(vare)* X1 
X2 ~~(vare)* X2 
X3 ~~(vare)* X3 
X4 ~~(vare)* X4 
X5 ~~(vare)* X5 
X6 ~~(vare)* X6 
'
Xlat<-lavaan::growth(modelXlat,data)
Xlat.fit<-round(fitmeasures(Xlat)[c("chisq","df","pvalue","cfi","srmr","rmsea","aic","bic")],digits=2)

遷移プロット

• plot.growth関数(自作)

source("script/plot.growth.R")
a<-plot.growth(fit.all, type="no")+theme_bw()
b<-plot.growth(fit.all, type="lin")+theme_bw()
c<-plot.growth(fit.all, type="quad")+theme_bw()
d<-plot.growth(fit.all, type="latent")+theme_bw()

適合度

fits<-round(fit.all$fit[
            c("chisq","df","pvalue","cfi",
              "srmr","rmsea","aic","bic"),],digits=2)
fits[,2]<-Xlat.fit
datatable(fits,option=list(dom='t'))

ポチるおまけ

推定結果の出力

parm<-parameterEstimates(fit.all$lavaan$quadratic)
parm[,5:10]<-round(parm[,5:10],digits=3)
datatable(parm[c(28:30,37:39),],
          extensions = 'Scroller', 
          options = list(dom= ' t',
                         deferRender = TRUE, 
                         scrollY = 200, 
                         scroller = TRUE))

推定平均値の遷移プロット

source("script/plot.growth.R")
g<-plot.growth(fit.all, type="quad")
g+ylim(3,4)+theme_bw()

　コード全体

model1 <-'

#切片因子の設定
i =~ 1*X1 + 1*X2 + 1*X3 + 
1*X4 + 1*X5 + 1*X6

#傾き因子の設定
s1 =~ 0*X1 + 1*X2 + 2*X3 +
3*X4 + 3*X5 + 3*X6

s2 =~ 0*X1 + 0*X2 + 0*X3 +
0*X4 + 1*X5 + 2*X6


#切片と傾きの分散
i ~~ i ; s1 ~~ s1 ; s2 ~~ s2; 

#因子間相関
i ~~ s1 + s2; s1 ~~ s2

#因子平均
i ~ 1 ; s1 ~ 1 ; s2 ~ 1

#誤差分散
X1 ~ 0; X2 ~ 0; X3 ~ 0
X4 ~ 0; X5 ~ 0; X6 ~ 0
'

切片因子の設定

• 因子負荷を1に固定

#lavaan code
i =~ 1*t1+1*t2+1*t3+1*t4+1*t5

前半の傾きの設定

前半の傾き(s1)の因子負荷を
区分時点以降同値に固定

#lavaan model code
i=~0*t1+1*t2+2*t3+3*t4+3*t5+3*t6

後半の傾きの設定

後半の傾き(s1)の因子負荷を
区分時点まで0に固定

#lavaan model code
i=~0*t1+0*t2+0*t3+0*t4+1*t5+2*t6

その他の設定

• 切片と傾きの分散を自由推定
• 因子間相関を自由推定
• 因子平均を自由推定
• 誤差分散を0に固定

#切片と傾きの分散
i ~~ i ; s1 ~~ s1 ; s2 ~~ s2
#因子間相関
i ~~ s1 + s2 ; s1 ~~ s2
#因子平均
i ~ 1 ; s1 ~ 1 ; s2 ~ 1
#誤差分散
bmi1 ~ 0; bmi2 ~ 0; bmi3 ~ 0
bmi4 ~ 0; bmi5 ~ 0; bmi6 ~ 0

推定の実行

library(lavaan)
model1.fit<-lavaan::growth(model1, data)

フィットしてる？

fit1.m<-round(fitMeasures(model1)[c("chisq","df","pvalue",
"cfi","srmr","rmsea")],digits=2)
fit1.m<-t(as.data.frame(fit1))
print(xtable(fit1.m),comment=F,type="html")

• piecewise

• quadratic

前半と後半の傾きに違いがあるか？

• 傾き因子の平均に等値制約を置いたモデルと比較

  #傾き因子平均が等値
  s1 ~ (a)*1 ; s2 ~ (a)*1

# model1.2fit = 等値制約のモデル
anova(model1.2fit, model1.fit)

## Chi Square Difference Test
## 
##             Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)
## model1.fit  12 26252 26336 324.30                              
## model1.2fit 13 26251 26330 325.46     1.1551       1     0.2825

もうちょい頑張る

• 前半を２次曲線モデルに

model2 <-'
#切片因子の設定
i =~ 1*X1 + 1*X2 + 1*X3 + 1*X4 + 1*X5 + 1*X6

#傾き因子の設定
s1 =~ 0*X1 + 1*X2 + 2*X3 +3*X4 + 3*X5 + 3*X6
s2 =~ 0*X1 + 1*X2 + 4*X3 +9*X4 + 9*X5 + 9*X6
s3 =~ 0*X1 + 0*X2 + 0*X3 +0*X4 + 1*X5 + 2*X6

#切片と傾きの分散
i ~~ i ; s1 ~~ s1 ; s2 ~~ s2; s3 ~~ s3; 

#因子間相関
i ~~ s1 + s2 + s3; s1 ~~ s2 + s3 ; s2 ~~ s3 ;

#因子平均
i ~ 1 ; s1 ~ 1 ; s2 ~ 1 ; s3 ~ 1

#誤差分散
X1 ~ 0; X2 ~ 0; X3 ~ 0; X4 ~ 0; X5 ~ 0; X6 ~ 0
'

　どやさ

-今回の適合度

model2.fit<-lavaan::growth(model2, data=data)
fit2.m<-round(fitMeasures(model2.fit)
[c("chisq","df","pvalue",
"cfi","srmr","rmsea")],digits=2)
fit2.m<-t(as.data.frame(fit2.m))
datatable(fit2.m,options=list(dom="t"))

• さっきのpiecewiseモデル

• 圧倒的なfit感!!

圧倒的なfit感!!

　期間ごとに切片が異なるモデル

• 二つ切片因子を設定

model2.2 <-'
i1 =~ 1*X1 + 1*X2 + 1*X3 + 0*X4 + 0*X5 + 0*X6
i2 =~ 0*X1 + 0*X2 + 0*X3 + 1*X4 + 1*X5 + 1*X6
s1 =~ 0*X1 + 1*X2 + 2*X3 + 3*X4 + 3*X5 + 3*X6
s2 =~ 0*X1 + 0*X2 + 0*X3 + 0*X4 + 1*X5 + 2*X6

i1 ~~ i1
i2 ~~ i2
s1 ~~ s1
s2 ~~ s2

i1 ~~ i2 + s1 + s2
i2 ~~ s1 + s2 
s1 ~~ s2
i1 ~ 1
i2 ~ 1
s1 ~ 1
s2 ~ 1

X1 ~ 0
X2 ~ 0
X3 ~ 0
X4 ~ 0
X5 ~ 0
X6 ~ 0 '

こんな感じでプロット

成長曲線モデルまとめ

• 遷移パターンは複数のモデルを比較検討
• 切片のみ、線形、二次曲線、latent basis
• RAMpath パッケージで瞬殺
• 区切りが明確な場合はpiecewiseも試してみよう

2つの時系列データ

• ポジ感情とポジ認知の経時データ

Garland, E. L., Geschwind, N., Peeters, F., & Wichers, M. (2015). Mindfulness training promotes upward spirals of positive affect and cognition: multilevel and autoregressive latent trajectory modeling analyses. Frontiers in psychology, 6.

有効なモデル

• parallel process growth model
• auto regressive latent trajectory model
• latent change score model

1. 各変数に最適な潜在成長曲線モデルを検討
2. 潜在変数(切片と傾き)間の相関を推定

ALTM

２つの時系列データ

交差遅延効果モデル

成長曲線モデル

交差遅延効果+潜在成長曲線=ALTM

自己回帰のパス

-変数の安定性

交差遅延のパス

• 変数間の因果関係

　切片と傾き (全体の遷移パターン)

さ

し

せ

そ

た